Making big steps in trajectories 



Norbert Th. Miiller* Margarita Korovina^ 

Abteilung Informatik, FB IV CICADA 
Universitat Trier, Germany University Manchester. UK 

mueller@uni-trier . de Margarita. KorovinaOmanchester .ac.uk 



We consider the solution of initial value problems within the context of hybrid systems and empha- 
sise the use of high precision approximations (in software for exact real arithmetic). We propose a 
novel algorithm for the computation of trajectories up to the area where discontinuous jumps appear, 
applicable for holomorphic flow functions. Examples with a prototypical implementation illustrate 
that the algorithm might provide results with higher precision than well-known ODE solvers at a 
similar computation time. 



1 Introduction 



The central idea underlying hybrid systems is that of a system of differential equations (initial value 
problems, IVP) enhanced with the ability to do discontinuous jumps, simular to finite automata. Unfor- 
tunately, from the viewpoint of TTE (e.g. Ill) discontinuity implies non-computability, which has been 
investigated in detail in [4], e.g.. Nevertheless, the importance of these systems forces us to deal with 
them and provide the best solutions possible. 

There do exist many software tools for hybrid systems (see lIH e.g.); however, almost all of them are 
based on double precision numbers; a notable exception is Ariadne lUl using generic programming, 
it is hence prepared for other implementations of real numbers. 

One basic aspect of the hybrid systems is their evolution in time, i.e. the computation of trajectories. 
In this paper we present an algorithm (implemented using the iRRAM package) for efficient and arbi- 
trarily precise solutions of the underlying IVPs up to the area where discontinuous jumps appear. There 
are two reasons for using high precision solutions: firstly, low precision might lead to incorrect assump- 
tions about the location of these jumps points; and secondly - perhaps unexpectedly - high intermediate 
precision can sometimes increase the efficiency. 

To illustrate the second aspect, consider the well-known Runge-Kutta methods used for the solution 
of IVPs. These methods are of fourth order, i.e., the error depends on a bound on higher derivatives of 
the solution as well as on the fourth power of the step width. Although they are a reasonable choice if 
applied to double precision numbers, they will not always be optimal for higher precision solutions. 
As the step width has to be chosen according to the desired precision of the solution, methods with a 
fixed order lead to the number of steps growing exponentially in the number of bits of the solution. If the 
order can be chosen dynamically and arbitrarily high, significantly fewer steps associated with a much 
bigger step width are possible, which can lead to a polynomial time complexity lITOll . 

Differential equations have been addressed numerous times in computable analysis, for example 
see ||2lil5|, where general questions of computability are addressed. A very important related result 
can be found in [7 |: The computation of solutions of differential equations is closely related to the 
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problem '#P=FP' from discrete complexity theory. Ttiis immediately implies that for general IVPs we 
cannot expect to find algorithms that run in polynomial time. For special IVPs, on the other side, it 
is well-known that the solutions are computable in polynomial time. A very detailed analysis of the 
resulting complexity for one-dimensional solutions can be found in ifTOl : If the flow function of the IVP 
is holomorphic and computable in polynomial time, then we are able to use methods with arbitrarily high 
order to solve the IVP and get a polynomial time solution. 

In this paper we will take the result from [10] and generalise it to IVPs of arbitrary finite dimension. 
This generalization then is used as the fundamental part of a new algorithm for the computation of 
trajectories in hybrid systems. In addition to a discussion of the theory behind our approach we actually 
present a prototypical implementation in the iRRAM software package |[TTl[T2ll . As an example, we use 
a rather simple linear differential equation where even an analytical solution is known. This allows us 
to compare our implementation with conventional IVP solvers, where our prototype has already shown 
unexpected efficiency at an always superior precision. 

Implementations of IVP solvers in exact real arithmetic, giving arbitrary precision results, are very 
rare. A prototypical implementation mentioned in [5 ] can hardly be useful in practice, as it seems to be 
based on the explicit construction of the solution using piecewise linear functions. This will necessarily 
lead to a complexity that is exponential in the precision of the solution. Vaguely similar approaches 
shown in the tutorial section during the CCA 2009 conference in Ljubljana were already unable to com- 
pute more than 4 decimals of the integral Jq f{t)dt for the simple function /(?) = t^. This leads us to 
assume that in this paper we actually present the first usable implementation for IVPs using exact real 
arithmetic. We should mention here that IVP solvers using interval arithmetic (hence also correct, but 
not arbitrarily precise) are well-known, for an overview see ||T?I . e.g.. 

2 Hybrid Systems 

A hybrid system can be defined as a tuple H = (2,X,D, G,F,R) consisting of a finite index set Q, 
a continuous state space X = VjqfzqXq, a collection of invariant domains D = {D^j^eg, Dq C Xq, a 
collection of guard sets G = {Gq}qeQ, Gq C Xq, a collection of continuous dynamics or flow conditions 
F = {Fq}q^Q defining differential equations in X, and a collection of discrete dynamics or reset relations 
R = {Rq qt^o], Rq qi ^ Xq X Xq' , scc e.g. [14J. The part that we are addressing in this paper are the 
flow conditions Fq that lead to trajectories in a component Xq of the state space. Whenever such a 
trajectory enters the guard set Gq, a discontinuous jump according to R might happen, quite similar to a 
non-deterministic automaton with state space Q. 

In the following, we will consider single trajectories of such hybrid systems. Our goal is not to 
discuss their computability or to formally consider their (theoretical) computational complexity, but we 
want to get a usable implementation that computes such trajectories within a component Xq until they 
enter the guard set Gq. As we do not attempt to use the reset relations, we will not really need Q, and in 
consequence we will omit all references to Q in the following. 

Below we will shortly describe the data structures we use; they will implicitly impose restrictions on 
the hybrid systems we can address: 

• State space X: We will use X = M x M'^' for an integer dimension (i G N. In the implementation, we 
use (dynamically sized) vectors of real numbers, d is then not given explicitly, but can be derived 
from the size of the vectors. The first component of X will always be used as the time parameter. 

• Invariant domain D: For the flow conditions F that we are able to address at the moment, it would 
be very artificial to use a proper subset of X here. For simplicity, we only use D = X. 
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In the following we will nevertheless mention where extensions for a non-trivial D would be nec- 
essary. In general, the distance function dx\D{^) from vectors G X to the exterior of D should 
be computable, implying that D should be computably open. 

• Flow conditions F: The basic property of the flow conditions is that there are d functions Fy : 
C(M X M'^) — )• M defining differential equations. We will assume that the Fy are holomorphic, 
hence they can be smoothly extended to complex arguments. In the cases we are able to deal with 
now, we will even use that the Fy are holomorhpic on the whole set C^^^ 

If the Fy are only holomorphic in a restricted area (with D C Dc), the distance function 
^C''+i\i»c ('= ) froni vectors G D to the exterior of Dc should be computable. 
At the moment, our implementation is restricted to flow functions Fy that are in fact polynomials 
(in d+l variables and with computable coefficients). We will use the following data to get access 
to the relevant properties of F: 

- Maximal degree n of the polynomials 

- Single coefficients of the polynomials, i.e., we have direct access to the d ■ {n + 1)^+' real 
numbers Cy^kj^ j^ for 1 < v < J and < ii,...,id,k < jj.. In applications, many of these 
numbers will be known to be zero. Our later examples of linear homogeneous IVPs will 
even be restricted to only d ■ {d + 1) values that may be non-zero. 

- A function Up on states {to,wo) € X and distances 5,£ S returning an upper bound for 
the flow functions Fy on the compact set 

C((fo,wo),5,£) := {{t,z) e C^+i : \t-to\ < 5 A |z-wo| < e} 

Such a function can obviously be computed directly using the coefficients, implying that Uf 
would be superfluous in a minimalistic setting. Nevertheless, we will later see that Uf plays 
a special role in the solutions and should also be helpful in cases that are not yet covered by 
our implementation. 

• Guard set G C X. G will be given by an algorithm to compute dci^), which implies that G is 
restricted to computably closed sets. Additionally, for a part of our algorithm we will use that its 
interior G" is computably open; so also the distance dx\G"{^) to the complement of G" must be 
computable. 

• Trajectories: Given an initial state ^ = {to,wo), our goal is to to determine how long the trajectory y 
through (fo, >vo) lives until it enters the guard set G; we want to find the first t >to such that y{t) is in 
G. Such a trajectory y is a vector (ji , . . . ,yd) of (possibly partial) real valued functions yy : CM — M 
on a real variable t that is usually interpreted as a 'time' parameter. As y touches {to,wo), y is (a part 
of) the solution of the IVP with flow conditions F and initial condition y{to) = wq. In our setting, 
it will be natural to extend this to complex variables t and to consider yy : CC — )■ C instead. 

3 Solving IVPs 

The basis of the approach we take is a well-known recursion that can be found in mathematical text- 
books like [3|. In [ 10| we used it to derive an algorithm showing polynomial time complexity for certain 
one-dimensional IVPs. We now generalise this result for higher dimensional systems. The polynomial 
complexity could also be provable in these cases, but in this paper we concentrate on an actual imple- 
mentation. 
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3.1 General solution for holomorphic flow functions 

Suppose we have d functions Fy : M x M'' — )■ M (l<v<<i) and a vector vvq = (wi^O; • • • ,Wdfi) such that our 
IVP has the form 



yy{t)=Fy{t,yi{t),...yd{t)) , 3'v(0)=Wv,o (1<V<J). 



(3.1) 



Note carefully that here we consider the special case where the initial condition {to,wo) is restricted 
to to = 0. The general case of specifying an arbitrary time to will be addressed later. 

The Picard-Lindelof theorem guarantees a unique solution on some interval containing to = if the 
function vector F and its partial derivatives are continuous on a region around {to,wo). This theorem 
can be used to get an iterative algorithm for the construction of solutions, like in ||5l. We will now use 
much stronger conditions: Our assumption is that the functions Fy are holomorphic, i.e., there exists a 
system of coefficients Cy^kJu-Jd such that for f € M near the origin, for x = (xi , . . .x^) G and for each 
coordinate V (l<v<J) we have 



Fv{t,x) 



<:,!l,...,rrfGN 



r r ■ • • • r'' • • y''' 
v,*:ji 1 • ■ ■ d 



(3.2) 



We assume that the circle of convergence is large enough to contain the initial condition (0,wo). (Cen- 
tering the circle of convergence at wq as well as using an arbitrary to will be addressed later.) Then 



A:,!l,...,(rfGN 



(3.3) 



As the Fy are holomorphic, the functions yv{t) are also holomorphic. This follows from 'local' versions 
of the Picard-Lindelof theorem in the complex plane, but we can also just use a Taylor series approach 
and prove that the radius of convergence is not zero (which is essentially done in section 4 below). So, 
let Qy^n be the corresponding sequences of coefficients, such that 



yv{t) = «v,n-f"- 

The coefficients ay o are already known from the initial condition, as we have 

Qv.o =3'v(0) =Wvfi . 



(3.4) 



(3.5) 



In the following, we address the coefficients a^^n with « > 0. From equation (3.3 1 we see that we need 
the powers {yy{t)y . The corresponding coefficients will be denoted by ay]„, so 



Comparing coefficients, we get the following recursion for the coefficients of the powers: 



(3.6) 



(i+i) 



n-th value of the sequence 1,0,0,0,. 
V (0 

ay J • ^y^n-y • 



(3.7) 



0<j<n 
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It is worth noting that each ay\ is determined by the values a^j with j < n. A reformulation of (3.3 1 
now leads to: 



yv{t) 



I 

k,ii,...,ide'M 

I 

k,ii,...,ije'M 



Cv.k. 



;i,...,irf 



III 

j(:,il,...,('rfGN£GN «i,«2,---.nd6N 
«l + ...+nrf+/:=<' 



feN A:,ni,n2,...,nrfeN!i,...,!rfGN 



Cv,i. 



Comparing coefficients with jv(0 = I (^+1) ' '^v.i+i • we see that 



k,ni,n2,...,n,ieN h,...,ij€N 
ni+...+nj+k=i 



Cv,k. 



Aid) 



(3.8) 



Remembering that we already have av,o = i^v,o and that the values a"!'^, needed for £ + 1 only depend on 
flj./i for jJ. <nj{<i), this is in fact a recursion in the coefficients. 

While for any given £ the outer sum (using ni+...+nd+k = ^) in (3.8 1 is finite, the inner sum (using 
/i, G N) nevertheless usually forces us to do an infinite summation. 

In the general case, trying to compute the coefficients this way would involve non-algebraic methods: 
To compute the sum in equation (3.8 1 we need to show that the coefficients Cy^k.i, y converge to quite 
fast, a similar situation to the summation of a power series. From the experience in 191, we make the 
conjecture that this is true and that a uniform polynomial time complexity of the coefficients should lead 
to polynomial time complexity of the sums l^c...... as well as of the sequence (a^). In this paper we do 

not want to generalise the lengthy proof for the one-dimensional systems given in [9J, but rather have a 
closer look at several special cases that lead to further efficiency. 



3.2 Zero vector as initial value 



If we not only use fo = but additionally restrict the initial value wq to (0, ... ,0), (|3.8[) can 



significantly. In this special situation we have = av,o = Wv,o = 0. In the recursion (3.7 1 this implies 



.(2) _ „(2) 



*v,l 



*v.O 



0, then further a 



(3) 
V.2 



,(3) _ .(3) 



'v,l 



'v,0 



etc., by induction we get av\ = for / > n. This 



DC simplified 



finally reduces (3.8 1 to the following form: 



1 



£+1 



k,ni ,n2,...,nj£N 0<r'i <ni,...,()<id<nd 
ni+...+nd+k=£ 



'V,k.i, 



J'd) ' 
d,nd 



(3.9) 



Equations (3.7i and (3.9 1 together provide us with a finite recursion scheme to compute all coefficients 
of the Taylor series. Using exact real arithmetic like the iRRAM software this can be implemented quite 
straightforward. 
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3.3 Non-zero initial time 

Additionally to an arbitrary wq we might also consider arbitrary time instants to for the initial condition 
wq = (wi.o, • • • , Wd,o) such that (to) = Wy.o- Usually, this is ignored as we are able to transform the flow 
functions Fy to 

Ev{t,xi,...,Xd) ■.= Fv{t + to,xi +wifi,...,Xd + Wdfi) 
and consider the following system: 

Zvit)=Ey{t,ziit),...Zdit)) , zv(0) = i\<v<d). (3.10) 

Let Zv{t) be the solution functions for this system. The original system using the Fy and yv{to) = Wv,o is 
then solved by the functions yv{t) defined as 

yv{t) := zv{t -to)+Wvfl ■ 

Unfortunately, we now need the coefficients of power series for the Ey, that is for Fy centered at {to,wo) 
and not at the origin (0, ... ,0). How to get these coefficients depends on how Fy is given. 

• If Fy is given via the series (centered at (0, . . . ,0)), we could do a re-expansion in the new center, 
which is essentially a composition of a power series and a linear function. Here 10 could be a good 
starting point, where the composition of two power series is considered (but just for one variable). 

• If Fy is given via an algorithm computing the function in a neighborhood of the initial value, we 
could try to develop Fy into a series, similar as in |9|. 

In any case, trying to compute the coefficients for Ey will force us to use non-algebraic methods and 
to evaluate quite complicated infinite sums. Our implementation is not yet mature enough to treat such 
complicated cases; in the following we look at cases where we can have arbitrary initial values without 
the need to do a complicated transformation of the power series. 

3.4 Autonomous linear systems 

A very important class of applications are linear systems of differential equations, i.e., each function Fy 
is actually linear in the arguments Xj : 

Fy{t,Xl,...Xd) =/v,o(0+/v,l(0 ■Xl + ■.. + fy,d{t) -Xd . 

In this case, the IVP is already reduced to a system ofri^ + n functions fyj : M — )• M. 

Note carefully that fyj do not need to be linear themselves. They still can be very complex: in 
our setting they would be holomorphic (but now in just one variable t). These linear systems lead to a 
coefficient system with 

Cv,Ui,-,id 7^ =^ ii + ... + id<l . 

Still infinitely many coefficients could be non-zero (varying with k), leading to the necessity of an infinite 
summation in ( |3.8[ ). 

In applications, a further property of the IVPs under consideration might be helpful: often the systems 
are autonomous, i.e., the flow conditions do not depend on the time parameter t. This implies that each 
function fyj actually has to be constant, so in this case our coefficient system satisfies the following 
restriction: 

Cv,kM,.j,^0 =^ ii + ... + id<lAk = . (3.11) 
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In this case we only have to consider (f- + d numbers (instead of functions) defining the IVP. 



Together with the recursion for the from equation (|3.7l), we get the following finite(!) recursion 



scheme for the computation of all coefficients of the solution of the IVP: 



1 



'+1 



«l ,n2,...,nj&i 
ni+...+n^=l 



I 

'l,'2v,'rfG{0,l} 
(l+...+(rf<l 



.. .-a 



d,nj 



(3.12) 



(0) 



This formula can be further reduced using that the values ay j, in (3.7 ) are very simple: As soon as a term 



in (3. 12 1 contains a component a^'^/ with ij = and nj > 0, the product is zero. On the other hand, the 
condition /i + . . . + < 1 implies tiiat at most one index ij can be non-zero. 
In consequence, for Z > we get the recursion 



1 



T I 



(1) 

Cv,0,0,0,...l,0,- -,0 ■'2^7 



(3.13) 



Additionally, for / = 0, we have one further term: 



'S'v.l — Cv,0,0,...,0+ 52 



Cv,o,o,o,...i,o,...,o'fl^j|d 



(3.14) 



Please note that for these autonomouse linear systems, the power series for the flow functions trivially 
have an infinite radius of convergence. Furthermore, the initial condition {Iq^wq) can now be arbitrary: 
we may use the transformation from the previous section, but only applied to the time parameter t (as we 
are able to treat non-zero wq directly). But then, due to the autonomous system, the coefficients for the 
transformed system and the original system are identical; we only have to solve 

Zv(t)=F^{t,zx(t),...Zd(t)) , Zv(0) = Wv,o il<v<d) . 

using ( 3.5|3.7|3.13|3.14 i. This gives us a power series for each Zv and we simply have use 

yvit) ■=Zvit-to) . (3.15) 



3.5 Nonlinear type, not autonomous, using multinomial flow functions 

From the previous section it is clear how to get a bigger class of IVPs that can be treated algebraically. 
Suppose there is some /i G N such that Cv^k,ii,...,ij is zero as soon as one of its indices is larger than 
jj.. Then the infinite sum 



in (3.8 1 reduces to just a finite sum ^ 



Additionally, 

re-centering the coefficients like in Subsection |3.3| is just a finite manipulation of polynomials. 

In consequence, the current version of our implementation contains an IVP solver based on the 
following summary of the considerations in this section: 



Suppose in (3.2 1 only a finite number of coefficients Cy 



are non-zero. 



Suppose that the initial value {to,wo) as well as the Cv,k,iu...,ij are uniformly computable. 
Then the power series av,n for the solutions of the IVP are uniformly computable. 
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4 Taylor Sequences and Bounds 

The previous section has shown how we can get access to the coefficients of the power series for the 
IVP solution. An important additional part of the evaluation of the IVPs is, of course, the necessary 
summation of these sequences. Here we can obviously not avoid to compute sums of infinitely many 
values. The computational complexity of such a summation has been addressed for example in f9l: if 
the coefficients of the series are uniformly computable in polynomial time, then the sum function also 
has polynomial complexity in the interior of the circle of convergence. To show this, a very detailed 
consideration of all intermediate rounding errors was necessary. This would also be necessary in an 
implementation if we use a traditional multi-precision package for the computation. Fortunately, an 
implementation using exact real arithmetic is much simpler in this regard, as the software package is 
able to deal with all these cumbersome details by itself and we can instead concentrate on the more 
important issue: the truncation error coming from using a finite summation instead of an infinite one. 

So consider a sequence (at) of Taylor coefficients together with a radius /? G M such that Y^a^ ■ 
converges absolutely (in the complex plane) for any z G C with \z\ < /? to a function /. Please note that 
essentially a simple linear transformation of the argument is sufficient if we have a series of the form 
'Lafiz-zof. 

In finite time, we are only able to compute partial sums L^^q '^k • z'^ from the sequence, leading to 
truncation errors of | cik- z'^l- However, we additionally just need a computable upper bound for 

this truncation error that converges to zero with increasing n to implement the infinite sum in exact real 
arithmetic. Suppose we have access to an upper bound M G M for |/(z)| on {z G C : |z| = R}. Using 
the Cauchy integral formula, we see that < M -R^" holds for any n. This gives rise to the following 
expUcit error formula, valid in case of |z| <R: 

OO CO OQ 

I I au-^\ < £ |«,|-|/| < I M./?-^.|z^| 

k=n+\ k=n+\ k=n+\ //I 1 \ 

I OO /I i\k Ayr D /I K^-'^) 

\z\\ v / \z\ \ M-R I \z\ 



So, for an approximation with a truncation error of < 2^ we just have to add all the terms a„ • z" until 

, n+\ 

becomes smaller than 2P. 



M R / \z 



R-\z\ \R 

The key to this summation is knowledge about pairs {R,M) with |/(z)| < M on {z G C : |z| = R}- 
In our case, the functions / under consideration are solutions yv of IVPs; we can construct such bounds 
using the underlying flow conditions F, if we have access to a bound for F. 

For any compact complex set C C C^^^ the maximum IJ.{F,C) := max^gc.i<v<rf l^v(i?)| exists as 
soon as the Fy are continuous on C; /X as a functional is even computable using standard representations 
for its arguments. For {to,wo) and arbitrary 5,e > we now consider the special complex neighborhood 

C := Cc{{to,wo),d,e) := {{t,z) G : \t -to\ < 5 A\z-wo\ < s} . (4.2) 

For the invariant domain D = Z we obviously have C C Z; for smaller D we would have to restrict 5 and 
e to 5,£ < dx\D{to,WQ) to ensure C C X. 

Unfortunately, using a general algorithm to compute fj. {F, C) from its arguments F and C could be 
very time consuming. Furthermore, we obviously do not need the exact maximum but only a (good) 
upper bound. So, for an efficient implementation, it is important that the flow conditions F are not only 
given with algorithms computing Fy or the coefficients Cv,k,ii,...,id, but also with an additional function 
Up v^ithin IJ.{F,C) < ((fo,vvo),5,£). 
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As long as a trajectory (considered as a function of a complex time variable!) does not leave this set 
C, it cannot change faster than given by the bound /i(F,C). So if we take the pair {R,M) given by 

U := UF{{to,wo),5,e) , 

R := mm{5,£/U} , (4.3) 
M := \wv\+R-U , 

we have \yv (z) | <M for all z G C with |z - | < /?. 



5 Meeting the Guard 

With a combination of the two previous sections, we are able to compute the solution of many inter- 
esting IVPs (at least on a small interval). Concerning the intended application to hybrid systems, we 
additionally want to find the point where a trajectory hits the guard for the first time. Using the function 
A{t) = dG{t,y{t)) of the distance between the guard set and the trajectory at time t, this is equivalent to 
finding the supremum to ■= sup{f > to \ (Vf',fo <t' < t)A{t') > 0}. 

To find to, we restrict ourselves to computably closed guard sets G, so the distance function dc is a 
computable real function with G = dQ^{{0}). In this case, tg is left computable Q. In the following we 
present an algorithm to actually approximate to from the left; later we also discuss an additional part of 
the same algorithm that can sometimes dehver approximations from the right, so that in the well-behaved 
case we even get a computable to- 



5.1 Steps 

Using the initial condition vvo at to we want to find a ti such that A{t) > for all t G {to,ti). First we have 



to choose 5,e > such that the (complex) neighborhood Cci{to,wo),5,£) from (4.2i lies in Dc. Then 



we compute U, R, and M as given in (4.3 1; let R' = R/2. 
Let ti :=to + si for 

si:=min{A{to)/U,R'} . (5.1) 

Using R and M, the previous sections allow us to compute the value vvi := y{ti) of the trajectory at ?i. 
Additionally we are sure that U bounds the flow functions in that region so that additionally between to 
and ?i the trajectory may not touch the guard; n.b. using R' = R/2 above is not crucial, we only need to 
be sure that si <R for (|4?T i. 



We may continue this process, having two options: we may extend the solution obtained so far 
using the already known Taylor series (algorithmically quite inexpensive, but only giving quite small 
extensions) or we may determine new Taylor series (algorithmically expensive, but leading to bigger 
extensions). These two options will be called 'small steps' and 'big steps'. 

1. 'small steps': If ti is still sufficiently smaller than to +R', we may determine the distance A{ti) 
between the guard and (fi,vvi) to find a new step size ^'2 = min{A{ti)/U,R' — si} and let t2 := 
ti +S2- We can compute '■= y{t2) and may even iterate this process leading to sequences {sj) of 
step sizes and (f,) of time instants with < to- 

Please note that here we are able to use the Taylor coefficients over and over again that have been 
computed from (?o,h'o). As always < it might be necessary to add further coefficients or 
to provide them with higher precision. But as we approach the guard, the step width — ti will 
converge to 0, so quite often we will need only a few (or even none) new additional coefficients. 
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2. 'big steps': As the resulting grow towards tc, the summation of the series (centered in to) gets 
increasingly difficult. But of course, we are free to use any (f,-, w,) as a new initial condition. This 
involves the necessity to compute new Taylor coefficients for y at the new center (f,-, w;). 

The computation of these coefficients is quite time consuming, but on the other hand, afterwards 



the evaluation is faster again because of an improved truncation error (4.1 1. Please note that now 
we will also have to recompute U, R, and M. 

The optimal point for the re-centering of the problem surely depends on the flow conditions F and 
also on the guard G. Later we will briefly mention a corresponding heuristic. 

Using this computation (in small or big steps) of points on the trajectory, we get an approximation of to 
from the left. 



5.2 Traversing the guard 

When additionally trying to approximate to from the right, we are in a situation similar to the computation 
of roots of functions, as tc = mm{t > to \ A{t) = 0}. The comprehensive analysis of root finding in O 
helps to identify conditions we should use in order to allow a successful computation of tc. In general, we 
will not be able to check whether these conditions are met, so we can only try to find right approximations 
to to- 

The most helpful condition for root finding is that the underlying function should 'change sign'. 
As the distance function do (and hence also A) we used so far is non-negative, we should have further 
information on the guard set G: If additionally the interior G" of the guard is computably open, the 
complement X \ G" is computably closed and the distance dx\G'' is computable, too. Then Yg{^) '■= 
dci^) — dx\G''{^) is zero only at the border of G. So instead of A(?) = dG{t,y{t)) we should rather use 
r{t) := YG{t,y{t)). This of course requires that the trajectory can be extended into the interior of G, which 
is trivial if the invariant domain is D = X. Otherwise, the hybrid system must be given accordingly. 

If r(?) really changes sign at ?g> the trajectory at ?g is not a tangent to the border of G. Unfortunately 
a corresponding test is uncomputable in general. At the moment we can only assume that the hybrid 
system and the initial condition are such that the sign change happens. If not, our algorithm will just 
compute the left approximation from the previous subsection. 

On the other hand, if r{t) changes sign at tc, then r{t) > for t < tG and there is an £ > with 
r{t) < for f e {tG,tG + £)■ Standard search methods can now be used to compute tG', here we want 
to propose a more elaborate method: A usual assumption in the world of double precision arithmetic 
is that the distance between the trajectory and the guard is differentiable with a derivative significantly 
different from zero at tG- In the following we translate this optimistic approach to the world of exact real 
arithmetic. 

So suppose the guard has a smooth surface, such that the distance to the guard is a continuously 
differentiable function (in R^+'). As the trajectory is differentiable, the function A{t), defined above as 
the distance between the guard and the trajectory at time t, will be differentiable, too. This allows us to 
estimate the time the trajectory still needs until traversing the guard: If and f,- are consecutive time 
instants where we evaluate the trajectory (outside of G), then di := (A(f,) — A ))/(?,• — is the 
derivative di = A.{Q for some time point between and f,-, due to the mean value theorem. Hence we 
may estimate tG ~ ti + p,- for p,- := A(f,)/5,-, and f,- := + 2 • p,- is a candidate where r(f,) < might be true. 
We only have to be careful when di is (almost) zero, for which we implemented the following algorithm: 
If our goal is to find tc with an error of at most 2^" , then we first check whether surely p, < 2^"^^; 
only in this case we really compute p, and r(f,). So, together with strictly monotonic increasing left 
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approximations for to, we can also construct a candidate list {tj) of possible right approximations, all 
satisfying tj < tj + 2^". 

If we know for sure that r{tj) < for such a candidate tj, we know tj < to < tj', so we have an 
approximation to to with error 2^" allowing us to terminate the algorithm. 

Unfortunately, we cannot check easily whether r(tj) < 0, as such a test is not computable. As a 
substitute for this un-feasible test we use multivalued tests whether r{te) < -1-'' or T{t() > -2'-*^ with 
a precision k. These tests can be computed in finite time and they are applied as follows: whenever we 
compute a new pair {tj,tj) we check all previously computed candidates t£ (i.e. for i < j) again, but now 
with higher precision k := j. As j goes to infinity, we will eventually find any f/ such that r(f^) < 0; as 
soon as the first is found, our algorithm stops with success. 

The efficiency of these repeated tests is greatly improved, if we remove all those candidates tj from 
the list where we found (in a similar way) that r{tj) > 0. Additionally, we may remove all candidates tj 
from the list, as soon as a value p; is again larger than 2^"^^ 



6 Prototypical implementation 

We used the iRRAM package to implement a prototype for the proposed algorithm, whose core structure 
uses dynamically constructed functions of types FUNCTION<int , vector<REAL> > (for sequences of 
real vectors) and FUNCTION<REAL, vector<REAL> > (for vector- valued functions on the real numbers). 
The implementation of such function objects in an imperative language like C++ has been described in 
lfT2l . it is based on a lazy evaluation technique. Thus we avoid the necessity to implement explicit (and 
computationally very expensive) representations for functions and sequences given in |,2jil5J, e.g.. 
Using corresponding constructors 

• a=ivp_solver_simple (w,F) yielding a vector power series a for flow conditions F and an 



initial condition w implementing equation (3.12 1 



• f =taylor_sum(a,R,M) yielding the (vector-valued) sum function f for a (vector-)power series a 
and corresponding radius R and bound M, and 



• w=f (bs) evaluating f at bs with bs<R as a limit using equation (4.1 1 to control the truncation 
error, 

the core of the implementation is essentially just the following loop of 'big steps' interspersed with 'small 
steps' as mentioned in the previous section: 

do -[ // big steps 

a = ivp_solver_simple (w,F); 

. . . compute R,M ... 
f = taylor_sum(a,R,M) ; 
do -[ // small steps 

. . . compute a step size from the distance to the guard . . . 
. . . accumulate the step size in a variable s ... 
. . . evaluate f (s) . . . 

. . . try whether a sufficiently good approximation has been found . . . 

... if yes : stop . . . 
} until ( s is large enough for a big step ) 
w= f(s) 

} 
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As the evaluation of f (s) is a core part here, it is important to get a reasonably efficient implemen- 
tation of the Taylor summation. The existing limit operators in the iRRAM package were not fitting, as 
they were not yet applicable to FUNCTION objects (they were essentially only usable for predefined algo- 
rithms). Additionally, the general heuristic of the iRRAM (that tries to compute Umits with the maximal 
used precision) lead to an enormous waste of time; a new limit operator based on < \4.l \ had to be added. 
All other necessary operations were already present in the published version of the package. 

As a first benchmark we used the simple linear system 

ji(0 =3^2(0 ; j2(0 = -3^1 (0+ 0-02 -^2(0 with fo = 0,wo = (0,l) . 

Without the term ciooi = 0.02, the solution y would simply be the pair (sin, cos) ; with the additional term 
we still have an oscillation, but with a growing amplitude. 

As guard set we chose G = { {t,xi,X2) G | xi < —2}. Here the question was simply to approximate 
the first tc where yiitc) = —2 (we found to ~ 73.5422061995...). The size si of the small steps was 
chosen as in (5.1 1; whenever the accumulated small steps grew larger than rmn{~^/s\ ■ \/^,R'}, a big step 



was made. This bound of the big steps was chosen heuristically as an attempt to match the much higher 
complexity of the IVP solution at big steps with the more frequent Taylor summations at small steps. 

The following graph shows an 3d-plot of the resulting trajectory constructed from the (linearly inter- 
polated) points {ti,Wi). 
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The following table shows a few results of computations with this IVR Its interpretation is as follows: 
To approximate to with an error of at most 2^", the software chose a working precision of 2^^, using b 
big steps (re-evaluations of the IVP), s small steps (evaluations of the Taylor sum) with a maximal index 
of Imax (working with an order of imax) and took time t (on an AMD Athlon 64X2 Dual Core Processor 
4600+). 
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0.297s 
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12 
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5.42s 


10000 


11787 


14 


20361 
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To compare our results we used the IVP solvers from the popular high-level language octave, that is pri- 
marily aiming at numerical computations (www.gnu.org/software/octave), in order to solve the above 
IVP. We applied them just to approximate the trajectory starting from = up to 73.543. Only between 
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73.542 and 73.543 we tried to find tlie point wliere it dropped below -2 (without even trying to verify 
tliat tliis was tlie first solution). Witliin a few milliseconds, a naive application of the solver gave a result 
near = 73.54225. As only 6 decimals were in common with our result of to = 73.5422061995..., we 
tried less naive ways, which initially produced the same result tQ. Being convinced from the correctness 
of our own implementation, we continued playing with the octave solver; with further variations of its 
parameters applied in a quite elaborate way, we were able to get results different from both to and t^. 
The best combination we found resulted in 73.542208, now with 7 correct decimals, but within 0.7^ of 
computation time. As the results from the octave solver quite erratically jumped around 73.5422 with 
further variations of the parameters, we believe that more than 6 decimals precision cannot reliably be 
expected. Our conclusion from these experiments is that solving IVPs might be an area where exact 
real arithmetic can actually compete with ordinary double precision arithmetic in terms of speed and 
precision. 

To illustrate the effect of varying distances |fG — fo| on our algorithm, we removed the perturbating 
coefficient 0.02. Additional we chose the guard set Gn = {{t,xi ^xj) € M'' | ? > Tj} for a given r\ and just 
printed 9 leading decimals of y\ [to^ )■ As Iq^ = 1], this setting transformed our algorithm into a slow (but 
still exact) method to compute sin(7]). The results in the following table show that further reductions in 
the error propagation in our software are necessary before it can really be applied for larger ranges of 7] . 
Again we compared our results with the octave IVP solver, which was significantly faster for larger 7] 
but had problems with its precision again. 
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7 Summary 

In this paper we analysed a recursive method for ODE solving with an emphasis on algorithmic applica- 
bility. As the method consists of two basic parts, the construction of Taylor series and their subsequent 
summation, we were able to adopt it to special requirements found in hybrid systems. 

Based on our benchmarks we conjecture that a closer analysis of the complexity of the algorithm 
will show that for given poly-time computable F and initial values the computed value to has complexity 
polynomial in the precision, if the prerequisites of the algorithm are given. An open question is, whether 
this very specific aspect of precision-oriented complexity can also be expressed in a uniform way de- 
pending on F or at least on the initial value. It might be much easier to construct a dependency on the 
actual value of the derivative of A at to- Additionally, the effect of the heuristic for the relation between 
the 'small steps' and the 'big steps' on the complexity is also worth studying. 

Apart from these questions of computational complexity, further detailed considerations for the cases 
of non-linear equations and esp. for the non-algebraic solutions are important, as well as for the influence 
of the dimension of the state space onto the efficiency of the algorithm. 

Considering the trajectories within a single component of the state space is obviously only an initial 
step into hybrid systems. The far goal is to improve the efficiency of reachability analyses: for trajectories 
starting in a given subset of states we want to know all reachable states, including those induced by the 
jumps between the different components of the state space. Of course, treating single trajectories like 
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we did it in this paper can be extended to sets, as computability implies effective continuity. Our last 
benchmark concerning the solution of an IVP on longer time intervals however shows that the modulus 
of continuity derivable from our algorithm does not yet allow an efficient set-oriented evaluation. In 
the near future, we want to consider the influence of Lipschitz properties in the dependency between 
the initial condition and the final state at to to improve this. Additionally, we will address efficient data 
structures for closed sets together with corresponding set- valued functions. 

References 

[1] A. Balluchi, A. Casagrande, P. Collins, A. Ferrari, T. Villa & A. L. Sangiovanni-Vincentelli (2006): Ari- 
adne: a Framework for Reachability Analysis of Hybrid Automata. Proceedings of the 17th International 
Symposium on Mathematical Theory of Networks and Systems Kyoto, Japan, 24-28 July 2006. 

[2] Vasco Brattka, Peter Hertling & Klaus Weihrauch (2008): A Tutorial on Computable Analysis. In: S. Barry 
Cooper, Benedikt Lowe & Andrea Sorbi, editors: New Computational Paradigms: Changing Conceptions of 
What is Computable, Springer, New York, pp. 425^91. 

[3] I.N. Bronstein & K.A. Semendjajew (1985): Taschenbuch der Mathematik. Verlag Nauka, Moskau, BSB 
B. G. Teubner Verlags-Gesellschaft, Leipzig. 

[4] P. Collins (2008): Semantics and computability of the evolution of hybrid systems. Technical Report, Centrum 
Wiskunde & Informatica, Amsterdam MAS-R080L 

[5] A. Edalat & D. Pattinson (2007): A Domain-Theoretic Account of Picard's Theorem. LMS Journal of Com- 
putation and Mathematics 10, pp. 83-1 18. 

[6] J. van der Hoeven (2008): Fast composition of numeric power series. Technical Report 2008-09, Universite 
Paris-Sud, Orsay, France. 

[7] Ker-I Ko (1983): On the computational complexity of ordinary differential equations. Inf. Control 5^(1-3), 
pp. 157-194. 

[8] Jan Lunze & Frangoise Lamnabhi-Lagarrigue (2009): Handbook of Hybrid Systems Control. Cambridge 
University Press, Cambridge 

[9] Norbert Th. Muller (1993): Polynomial Time Computation of Taylor series. In: Proceedings of the 22th JAIIO 
-PaneJ'93, Part 2, pp. 259-281. Available at |http : //www . uni- trier . de/~mueller| Buenos Aires, 1993. 

[10] Norbert Th. Muller & Bernd Moiske (1993): Solving initial value problems in polynomial time. In: Proc. 22 
JAIIO - PANEL '93, Part 2, pp. 283-293. Available at , http://www.uni-trier.de/~mueller, Buenos 
Aires, 1993. 

[11] Norbert Th. Muller (2001): The iRRAM: Exact Arithmetic in C++. Lecture notes in computer science 2991, 
pp. 222-252. 

[12] N. Th. Muller (2009): Enhancing imperative exact real arithmetic with functions and logic. Technical Report, 
software presentation at the CCA 2009 conference , Ljubljana. Available at|http://www. uni-trier.de/ 
"mueller 

[13] NediaUco S. NediaUcov & Kenneth R. Jackson & George F. Corliss (1999): Validated solutions of initial value 
problems for ordinary differential equations. Applied Mathematics and Computation, 105(1) pp. :21-68. 

[14] A. J. van der Schaft & J. M. Schumacher (1999): Introduction to Hybrid Dynamical Systems. Springer- Verlag, 
London, UK. 

[15] Klaus Weihrauch (2000): Computable analysis: An introduction. Springer- Verlag New York, Inc. 



